Modeling-effort
Consultation notes with Daniela
Sarka Kadava
2024-10-17
Introduction
The experimental design is described in the preregistration here: https://osf.io/m9xep (@Daniela, you can skip the section Lab Equipment)
In this markdown, we will be modelling the causal relation between effort and correction to confirm/reject our hypothesis.
These are:
H1: In corrections, people will enhance some effort-related kinematic and/or acoustic features of their behaviour
H2: The enhancement will depend on similarity of the guesser’s answer and the original meaning. More similar answer will require/result in smaller enhancement (but still enhancement) than less similar answer.
Setting up the environment
library(here)## here() starts at E:/FLESH_ContinuousBodilyEffortlibrary(dagitty) # for dag
library(dplyr) # for data-wrangling##
## Attache Paket: 'dplyr'## Die folgenden Objekte sind maskiert von 'package:stats':
##
## filter, lag## Die folgenden Objekte sind maskiert von 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(lme4) # for linear mixed-effects models## Lade nötiges Paket: Matrixlibrary(tidyr) # for reshaping data (if needed)##
## Attache Paket: 'tidyr'## Die folgenden Objekte sind maskiert von 'package:Matrix':
##
## expand, pack, unpacklibrary(ggplot2)DAG
To make sure we know what variables should be in the model and why, we draw a DAG that illustrates our assumptions
Our relationship of interest is the causal effect of communicative attempt on effort.
Other assumptions include:
Personality traits (measured with Big5) will influence effort (e.g., people are more willing to put more effort if they are open-minded) and also communicative attempt (e.g., more extroverted people are better in this game, therefore they need less attempts)
Familiarity with guessing partner influences effort (ditto) as well as communicative attempt (e.g., friends are better in this game than strangers)
Similarly, participant will also directly infleunce effort and commAtt, because personality traits might not be capturing all variation
Expressibility of concepts is going to influence effort (e.g., more expressible concepts allow more enhancement - but could be also other direction) as well as CommAtt (e.g., more expressible concepts are more readily guessed and don’t need more attempts)
Similarly, concept will also directly influence effort and commAtt, because expressibility might not be capturing all variation
Modality (uni vs. multi) will influence the value of effort. We assume that in unimodal condition, the feature does not need to account for synergic relations with the other modality, and may carry the whole signal. In multimodal condition, however, there may be synergic relations that moderate the full expressiveness of this feature
Lastly, trial number (characterising how for one is in the experiment) could be hinting on learning processes through out the experiment, or - the other direction - on increasing fatigue
dag <- dagitty('dag {
Big5 [adjusted,pos="-0.823,0.657"]
CommAtt [exposure,pos="-1.033,0.028"]
Conc [adjusted,pos="-1.136,-0.848"]
Eff [outcome,pos="-0.102,0.025"]
Expr [adjusted,pos="-0.758,-0.850"]
Fam [adjusted,pos="-0.379,0.663"]
Mod_cat [adjusted,pos="-0.374,-0.850"]
Pcn [adjusted,pos="-0.589,1.214"]
TrNum [adjusted,pos="-1.686,-0.859"]
Big5 -> CommAtt
Big5 -> Eff
CommAtt -> Eff
Conc -> Expr
Expr -> CommAtt
Expr -> Eff
Fam -> CommAtt
Fam -> Eff
Mod_bin -> Eff
Mod_cat -> Expr
Pcn -> Big5
Pcn -> CommAtt
Pcn -> Eff
Pcn -> Fam
TrNum -> CommAtt
TrNum -> Eff
Conc -> CommAtt
Conc -> Eff
}')
plot(dag)## Plot coordinates for graph not supplied! Generating coordinates, see ?coordinates for how to set your own.
Synthetic data
We will now create synth data that will copy the relations we assume in our DAG. Assigning concrete coefficients will also help us to test our model - we should find exactly those causalities that we had in mind when creating these data
# Set seed for reproducibility
set.seed(0209)
# Define participants, total unique concepts, and modalities
n_participants <- 120
n_total_concepts <- 21 # Total unique concepts
n_concepts_per_participant <- 21 # Each participant works with 21 concepts
n_modalities <- 3 # gesture, vocal, combined
# Generate participant IDs
participants <- 1:n_participants
# Simulate Big5 personality traits (standardized between 0 and 1) and Familiarity (between 0 and 1) for participants
Big5 <- runif(n_participants, min = 0, max = 1) # Continuous values between 0 and 1
Familiarity <- runif(n_participants, min = 0, max = 1) # Continuous values between 0 and 1
# Create a matrix to hold expressibility values for each concept in each modality
expressibility_matrix <- matrix(runif(n_total_concepts * n_modalities, min = 0, max = 1), nrow = n_total_concepts, ncol = n_modalities)
# Randomly sample 21 unique concepts for each participant
final_data <- data.frame()
# Define a function to assign CommAtt and Eff for a single participant
simulate_participant <- function(participant_id) {
# Randomly sample 21 unique concepts from the total pool of 84
selected_concepts <- sample(1:n_total_concepts, n_concepts_per_participant)
participant_data <- data.frame()
trial_number <- 1 # Initialize trial number
for (concept_id in selected_concepts) {
# Randomly determine the modality for the concept
modality <- sample(c("gesture", "vocal", "combined"), 1)
# Calculate expressibility based on modality
expressibility_score <- ifelse(modality == "vocal", expressibility_matrix[concept_id, 1] * 0.6,
ifelse(modality == "gesture", expressibility_matrix[concept_id, 2],
expressibility_matrix[concept_id, 3] * 1.2))
# Determine Communicative Attempts based solely on expressibility, familiarity, and Big5
base_prob <- c(0.33, 0.33, 0.33) # Equal chance for 1, 2, or 3 attempts
# Modify probabilities based on familiarity, Big5, and expressibility
adjusted_prob <- base_prob * c(1 - Familiarity[participant_id], # 3 times for each
1 - Familiarity[participant_id],
1 - Familiarity[participant_id]) *
c(1 - Big5[participant_id],
1 - Big5[participant_id],
1 - Big5[participant_id]) *
c(1 - expressibility_score,
1 - expressibility_score,
1 - expressibility_score)
# Normalize the adjusted probabilities
adjusted_prob <- adjusted_prob / sum(adjusted_prob)
# Sample the number of communicative attempts based on adjusted probabilities
n_attempts <- sample(1:3, 1, prob = adjusted_prob)
# Loop through the number of attempts and increment CommAtt correctly
for (attempt in 1:n_attempts) {
# Calculate Eff for the first attempt
if (attempt == 1) {
Eff <- 1.15 * Big5[participant_id] +
1.10 * Familiarity[participant_id] +
1.20 * expressibility_score +
rnorm(1, mean = 1, sd = 0.5)
# Adjust Eff based on modality
if (modality == "combined") {
Eff <- Eff * 0.7 # Slight moderation for combined modality
}
}
# Adjust Eff for subsequent attempts
if (attempt == 2) {
Eff <- 1.15 * Big5[participant_id] +
1.10 * Familiarity[participant_id] +
1.20 * expressibility_score +
rnorm(1, mean = 1, sd = 0.5)
Eff <- Eff * 1.50 # Multiply effort by 1.50 for the second attempt
} else if (attempt == 3) {
Eff <- 1.15 * Big5[participant_id] +
1.10 * Familiarity[participant_id] +
1.20 * expressibility_score +
rnorm(1, mean = 1, sd = 0.5)
Eff <- Eff * 0.70 # Multiply effort by 0.70 for the third attempt
}
# Create row for each attempt
participant_data <- rbind(participant_data, data.frame(
Participant = participant_id,
Concept = concept_id,
Modality = modality,
Big5 = Big5[participant_id],
Familiarity = Familiarity[participant_id],
Expressibility = expressibility_score,
CommAtt = attempt, # Correctly set the attempt number
Eff = Eff,
TrialNumber = trial_number # Set trial number for this attempt
))
# Increment the trial number after each attempt
trial_number <- trial_number + 1
}
}
return(participant_data)
}
# Simulate data for all participants
for (i in participants) {
final_data <- rbind(final_data, simulate_participant(i))
}
# Preview the first few rows of the final data
head(final_data)## Participant Concept Modality Big5 Familiarity Expressibility CommAtt
## 1 1 3 combined 0.9599493 0.9219529 0.07567109 1
## 2 1 3 combined 0.9599493 0.9219529 0.07567109 2
## 3 1 3 combined 0.9599493 0.9219529 0.07567109 3
## 4 1 9 vocal 0.9599493 0.9219529 0.37757391 1
## 5 1 9 vocal 0.9599493 0.9219529 0.37757391 2
## 6 1 21 gesture 0.9599493 0.9219529 0.56918606 1
## Eff TrialNumber
## 1 2.624256 1
## 2 2.715052 2
## 3 2.650435 3
## 4 3.498290 4
## 5 5.902036 5
## 6 3.246798 6So now we have synthetic data where we exactly defined (using coefficients) what is the relationship between certain variables
- CommAtt -> Eff
The effort for second attempt is multiplied by 1.25 (Beta = 1.25) The effort for third attempts by 0.9 (ie decreases, Beta = 0.9)
- Fam -> Eff & CommAtt
Beta = 1.10 for Eff Negative effect on CommAtt
- Big5 -> Eff & CommAtt
Beta = 1.15 for Eff Negative effect on CommAtt
- Expr -> Eff & CommAtt
Beta = 1.20 for Eff Negative effect on CommAtt
- TrNum -> Eff & CommAtt
Beta = could be both neg (fatigue) or positive (learning)
Exploring structure (DP)
nrow(final_data)## [1] 5055final_data |>
janitor::tabyl(Participant, Concept)## Participant 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## 1 3 2 3 3 3 1 1 1 2 1 2 2 3 1 1 2 1 1 3 2 2
## 2 2 1 2 3 3 3 2 1 2 1 1 2 1 1 3 3 2 2 2 3 3
## 3 3 3 3 2 1 3 2 1 3 1 2 3 1 3 3 3 2 3 1 1 3
## 4 3 3 3 1 2 3 2 2 1 3 1 3 1 1 2 3 1 2 2 2 2
## 5 2 1 2 1 2 3 3 1 2 1 3 2 1 3 3 2 1 2 2 2 3
## 6 2 2 3 3 3 3 1 2 3 2 1 3 2 1 2 1 2 3 2 1 3
## 7 2 3 3 2 2 3 3 3 3 1 3 3 2 3 3 2 2 2 3 3 3
## 8 1 2 2 1 2 2 1 3 2 2 2 3 1 1 3 3 2 3 1 3 1
## 9 1 1 1 2 3 1 3 2 1 1 2 1 2 3 3 2 1 3 3 3 1
## 10 2 1 3 2 3 2 1 1 2 2 2 3 3 1 1 3 3 3 1 3 2
## 11 3 2 3 3 3 3 1 2 3 3 1 1 3 2 3 2 2 3 3 1 2
## 12 2 3 3 2 3 2 3 3 2 1 3 1 1 3 3 1 2 1 1 1 3
## 13 1 2 2 1 3 1 3 1 2 2 3 3 3 1 3 3 2 1 2 1 1
## 14 3 2 2 2 1 2 3 2 2 3 1 2 3 2 1 1 3 3 2 1 2
## 15 3 1 2 1 2 2 2 3 1 3 2 1 1 2 2 3 3 1 3 1 3
## 16 2 3 1 2 2 3 1 3 3 3 2 3 2 2 3 2 3 1 2 2 3
## 17 3 2 3 1 3 2 3 3 3 2 2 3 3 1 3 2 3 2 1 1 1
## 18 3 3 2 1 2 3 2 1 1 1 1 1 3 2 1 1 3 3 2 2 2
## 19 3 1 1 1 1 2 2 3 3 1 2 1 2 2 1 1 2 2 3 3 3
## 20 1 1 3 1 1 1 2 2 3 3 3 2 3 2 2 3 3 3 2 3 3
## 21 3 1 1 3 2 3 1 1 3 3 2 3 2 3 1 1 2 1 3 1 2
## 22 1 3 1 2 1 3 3 1 3 2 3 1 1 1 1 3 1 1 2 1 3
## 23 3 3 2 3 2 1 1 1 2 1 2 3 2 3 2 3 3 1 2 1 1
## 24 2 1 3 1 2 2 3 1 3 3 1 3 3 2 2 1 1 3 2 3 2
## 25 2 3 2 3 3 3 2 2 1 3 3 3 1 3 2 3 3 2 1 3 2
## 26 3 1 3 3 2 3 3 1 2 1 2 3 1 1 2 1 3 3 1 2 1
## 27 1 1 3 1 3 3 1 2 3 2 3 3 1 2 1 2 3 2 3 1 3
## 28 2 2 1 1 2 3 1 1 2 3 3 1 3 2 2 3 1 1 2 1 1
## 29 2 3 2 1 1 3 3 1 3 3 1 2 1 2 2 2 3 1 2 1 1
## 30 2 3 2 2 3 2 3 1 1 3 2 1 2 2 2 2 2 1 1 3 2
## 31 3 3 2 3 3 2 1 1 2 3 1 1 1 3 3 3 2 1 3 3 2
## 32 1 3 3 3 2 1 2 3 2 3 3 2 1 2 3 1 2 1 3 3 1
## 33 2 2 1 1 1 1 2 2 3 3 2 2 3 2 1 3 2 3 1 3 2
## 34 3 1 3 1 3 3 2 3 1 2 3 3 2 3 1 1 3 1 3 2 1
## 35 2 3 1 3 3 2 2 1 1 1 2 2 1 1 3 1 2 2 2 3 3
## 36 3 2 1 2 2 2 3 2 2 1 2 2 1 3 3 3 3 3 1 3 3
## 37 2 1 2 2 2 1 3 1 2 1 1 1 3 3 2 2 1 3 3 1 3
## 38 2 2 3 1 3 2 1 2 1 2 3 1 1 2 1 1 3 1 2 3 3
## 39 1 2 3 1 2 1 1 2 1 3 1 1 1 2 1 1 1 2 2 3 3
## 40 1 1 3 2 3 2 1 3 1 2 1 3 3 3 3 3 3 2 3 1 2
## 41 2 2 1 1 1 2 2 3 2 3 1 3 1 2 3 3 1 3 3 2 2
## 42 3 2 1 3 3 1 3 2 2 3 1 1 2 1 1 3 1 1 1 2 3
## 43 2 3 2 3 3 2 3 2 2 2 3 1 1 2 2 3 2 1 1 3 2
## 44 3 1 1 1 1 1 3 1 2 1 3 3 1 2 2 1 3 3 3 2 3
## 45 2 1 1 3 2 1 1 1 3 2 2 3 1 3 1 2 1 2 2 1 1
## 46 2 3 3 3 2 2 2 1 2 3 2 3 3 1 2 1 3 2 3 1 1
## 47 1 3 1 2 1 3 1 2 3 3 3 3 2 3 3 2 3 3 3 2 2
## 48 1 1 1 3 1 1 2 1 3 2 1 2 3 3 1 1 3 3 3 2 2
## 49 1 2 3 2 3 3 1 2 2 2 1 1 2 2 2 3 2 3 2 3 2
## 50 1 1 1 2 2 1 1 1 1 2 2 1 1 1 1 3 1 2 3 2 3
## 51 1 2 3 2 3 1 3 3 2 3 1 2 2 2 3 1 2 2 1 2 1
## 52 3 1 3 3 2 2 2 3 3 2 2 2 2 2 3 1 2 3 3 1 2
## 53 1 1 2 3 2 1 2 3 2 2 3 1 1 3 3 3 1 1 1 1 2
## 54 2 3 3 2 1 1 2 1 2 1 2 1 2 1 1 1 1 3 2 3 3
## 55 1 2 2 3 1 2 1 2 1 3 1 2 3 3 1 2 1 3 2 3 3
## 56 3 1 2 2 1 1 2 3 1 1 3 1 2 2 1 3 3 1 2 2 2
## 57 3 3 1 1 3 2 3 3 1 2 1 1 3 1 3 3 3 1 1 2 3
## 58 3 2 3 1 1 3 3 1 1 2 3 3 1 2 1 1 2 2 3 3 1
## 59 2 1 2 1 1 2 2 3 3 2 1 2 3 1 3 2 3 1 2 2 1
## 60 2 3 2 2 3 2 3 1 2 2 3 1 3 3 1 1 1 1 1 1 1
## 61 1 1 2 2 1 2 2 1 3 3 1 1 3 2 1 1 3 3 1 1 2
## 62 1 1 2 1 1 2 1 3 1 1 2 1 3 3 1 2 2 3 2 3 3
## 63 2 3 3 2 2 2 1 3 1 2 1 3 1 3 3 3 1 1 2 2 2
## 64 1 1 3 1 1 2 2 3 1 2 1 2 3 2 3 3 1 2 1 1 2
## 65 3 2 3 1 2 1 3 2 1 1 2 3 3 1 1 2 1 2 3 1 2
## 66 1 1 3 2 3 2 1 3 2 3 3 2 3 3 1 1 3 3 2 2 3
## 67 3 3 1 3 3 3 1 3 3 1 2 1 2 3 2 3 1 1 2 3 3
## 68 3 2 2 2 2 3 2 1 1 1 3 1 1 1 1 2 1 3 1 1 3
## 69 1 1 1 3 2 1 1 3 2 3 3 3 1 1 3 1 1 3 3 3 3
## 70 3 1 1 2 2 2 3 3 2 1 2 1 1 2 3 2 3 3 2 2 3
## 71 2 2 2 3 2 3 3 2 1 1 3 1 3 1 2 2 2 1 3 3 3
## 72 1 1 2 2 2 2 2 2 2 3 2 3 3 2 2 2 3 3 3 2 1
## 73 3 3 2 3 3 3 2 1 3 1 2 3 3 3 2 1 1 2 3 2 3
## 74 1 2 2 1 2 1 2 1 2 1 2 2 2 3 2 2 1 2 3 3 1
## 75 3 3 1 1 1 1 3 2 3 3 1 1 2 3 1 2 1 3 2 3 1
## 76 2 2 2 2 3 1 3 2 1 1 1 3 2 2 1 1 3 1 2 3 1
## 77 2 2 3 2 1 2 3 2 1 3 1 2 3 1 1 2 3 3 2 3 2
## 78 3 3 1 2 3 1 3 3 1 2 3 3 1 2 2 2 1 3 1 1 3
## 79 1 3 1 2 3 3 3 2 1 1 2 3 2 1 2 2 3 3 1 2 2
## 80 2 1 2 2 2 2 1 2 3 1 3 1 1 2 1 1 2 3 1 2 2
## 81 2 3 1 3 2 2 1 3 1 2 1 2 3 3 1 3 3 3 1 2 3
## 82 3 3 3 2 2 1 2 1 1 2 1 3 1 1 1 1 1 1 1 3 1
## 83 1 2 1 1 1 1 2 2 2 3 1 3 3 1 3 2 3 1 1 3 1
## 84 2 3 2 2 1 1 1 3 1 2 1 1 3 3 3 2 1 1 2 1 3
## 85 3 2 3 3 2 1 3 3 1 2 3 1 1 1 3 1 1 3 1 3 3
## 86 3 1 2 3 3 2 1 3 3 1 2 3 1 2 3 3 1 1 3 1 2
## 87 1 3 3 1 2 3 3 3 2 3 1 1 3 1 2 1 3 3 2 3 1
## 88 1 2 2 2 1 3 2 3 1 1 3 2 2 3 1 2 2 2 1 2 3
## 89 3 3 1 3 2 3 1 1 3 3 2 2 3 3 1 2 2 1 3 2 3
## 90 3 3 2 2 1 3 1 3 2 3 1 1 1 1 3 2 1 1 2 3 1
## 91 3 2 1 1 1 3 2 3 1 1 1 3 2 1 2 3 2 2 2 3 3
## 92 1 3 1 2 1 2 2 2 3 1 3 1 1 3 3 2 2 1 1 2 2
## 93 2 2 1 1 3 3 1 1 2 3 2 2 1 3 1 1 2 2 3 2 1
## 94 2 2 1 2 3 1 1 1 2 3 1 2 1 1 3 1 2 1 1 2 2
## 95 2 1 3 2 3 2 3 2 2 2 3 2 1 1 2 1 3 2 2 3 3
## 96 2 3 3 3 3 1 3 3 1 2 3 1 3 1 3 2 1 3 2 1 3
## 97 1 3 3 3 2 3 1 2 1 2 2 3 2 3 2 1 3 3 2 3 3
## 98 3 3 2 2 3 3 2 2 3 2 2 1 2 1 3 2 1 2 3 3 2
## 99 1 3 2 3 1 3 3 3 1 1 3 1 3 2 1 1 2 1 3 3 1
## 100 2 2 2 1 3 1 2 3 2 1 3 2 1 2 1 3 3 2 2 1 1
## 101 1 3 2 3 3 2 1 3 1 1 1 1 1 2 3 3 2 3 1 3 2
## 102 1 2 1 2 3 2 1 2 3 3 1 3 3 1 3 3 2 3 3 3 2
## 103 1 1 3 1 2 1 2 3 1 2 2 1 2 1 2 2 3 3 2 3 3
## 104 3 3 1 3 3 1 3 3 1 2 1 1 3 3 3 1 3 1 3 1 2
## 105 1 3 2 1 2 2 2 3 2 1 1 2 1 2 1 1 3 2 3 1 2
## 106 3 3 1 2 2 2 1 1 2 2 1 1 2 3 2 2 2 1 1 1 1
## 107 3 3 2 2 2 2 3 2 2 2 3 2 2 1 1 1 2 1 2 1 1
## 108 2 1 2 2 1 2 1 1 1 3 1 3 3 1 1 2 1 2 3 3 3
## 109 1 2 2 1 3 1 2 1 3 1 2 3 1 2 3 1 3 1 3 1 3
## 110 2 3 1 2 3 3 3 2 2 1 1 2 1 1 2 2 1 2 3 1 1
## 111 3 2 2 3 1 2 1 1 3 2 2 1 1 2 1 2 1 2 2 1 3
## 112 3 1 1 1 1 2 1 3 1 2 2 1 2 1 1 2 1 2 3 1 2
## 113 2 2 3 2 1 3 3 1 3 1 1 3 2 3 2 1 2 1 3 1 2
## 114 2 3 1 1 3 3 1 1 3 1 3 2 3 2 1 3 3 2 1 1 1
## 115 3 1 3 1 2 3 3 1 1 2 3 1 3 2 1 1 3 1 1 2 2
## 116 3 3 1 1 1 2 2 1 3 2 1 1 1 1 1 2 3 3 3 2 1
## 117 1 2 3 2 3 1 1 3 3 1 2 3 3 2 1 3 3 2 2 2 2
## 118 2 2 2 3 3 2 1 2 1 2 1 3 2 2 3 2 2 1 1 2 2
## 119 3 1 1 3 3 1 1 1 2 3 3 3 3 3 3 2 1 3 3 2 3
## 120 1 2 1 3 1 3 3 3 1 1 1 2 3 1 3 2 1 3 3 1 3Check visual
# Sample data (replace this with your actual data)
# df <- data.frame(CommAtt = ..., Eff = ...)
# Create a boxplot comparing Effort across different Communicative Attempts
ggplot(final_data, aes(x = as.factor(CommAtt), y = Eff)) +
geom_boxplot(aes(fill = as.factor(CommAtt))) + # Boxplot with fill based on CommAtt
labs(title = "Comparison of Effort Across Communicative Attempts",
x = "Communicative Attempts",
y = "Effort",
fill = "CommAtt") +
theme_minimal() +
theme(legend.position = "none") # Optional: remove the legend
Modeling
Now let’s model our synth data
Contrast coding
# Convert necessary columns to factors
final_data$CommAtt <- as.factor(final_data$CommAtt)
final_data$Modality <- as.factor(final_data$Modality)
final_data$Participant <- as.factor(final_data$Participant)
final_data$Concept <- as.factor(final_data$Concept)
final_data$TrialNumber <- as.numeric(final_data$TrialNumber) # Ensure TrialNumber is numericDP: Check contrasts of factors
contrasts(final_data$CommAtt) <- MASS::contr.sdif(3)
contrasts(final_data$Modality) <- contr.sum(3)/2Center trial number
final_data$TrialNumber_c <- final_data$TrialNumber - median(range(final_data$TrialNumber))
range(final_data$TrialNumber_c)## [1] -26.5 26.5range(final_data$TrialNumber)## [1] 1 54Standardize (z-score) all continuous predictors.
But if the measures of familiarity, expressibility, and big5 are on a rating scale we can just subtract the median again from these to centre them (and we don’t need to standardize).
final_data$Familiarity <- final_data$Familiarity - median(range(final_data$Familiarity))
final_data$Big5 <- final_data$Big5 - median(range(final_data$Big5))Z-score expressibility (because it’s continuous)
final_data <-
final_data |>
mutate(Expressibility_z = (Expressibility - mean(Expressibility))/ sd(final_data$Expressibility, na.rm = T))Simple linear mixed model
Before we do full Bayesian, we do quick LMM to see whether we are on good track
# Fit the linear mixed-effects model
lmer_model <- lmer(Eff ~ CommAtt + Familiarity + Big5 + Expressibility + TrialNumber_c + Modality + (1 | Participant) + (1 | Concept), data = final_data)
# Summary of the model
summary(model)
# Check model diagnostics
# Plot residuals vs fitted values
plot(model)
# Extract coefficients
coefficients <- summary(model)$coefficients
print(coefficients)
# If you want to save the model output
saveRDS(model, here("06_Modeling", "models", "linear_mixed_effects_model.rds"))lmer_model <- readRDS(here("06_Modeling", "models", "linear_mixed_effects_model.rds"))
summary(lmer_model)## Linear mixed model fit by REML ['lmerMod']
## Formula: Eff ~ CommAtt + Familiarity + Big5 + Expressibility + TrialNumber +
## Modality + (1 | Participant) + (1 | Concept)
## Data: final_data
##
## REML criterion at convergence: 8420.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6533 -0.6248 -0.0245 0.6415 4.9419
##
## Random effects:
## Groups Name Variance Std.Dev.
## Participant (Intercept) 0.013366 0.11561
## Concept (Intercept) 0.001346 0.03668
## Residual 0.297073 0.54504
## Number of obs: 5067, groups: Participant, 120; Concept, 21
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.0748620 0.0449111 1.667
## CommAtt2 1.1706126 0.0171346 68.319
## CommAtt3 0.1163042 0.0217732 5.342
## Familiarity 1.1668413 0.0483836 24.116
## Big5 1.2090271 0.0460398 26.260
## Expressibility 1.1281766 0.0294578 38.298
## TrialNumber -0.0011358 0.0006236 -1.821
## Modalitygesture 0.9456262 0.0203315 46.511
## Modalityvocal 0.9134778 0.0206897 44.151
##
## Correlation of Fixed Effects:
## (Intr) CmmAt2 CmmAt3 Fmlrty Big5 Exprss TrlNmb Mdltyg
## CommAtt2 -0.148
## CommAtt3 -0.104 0.320
## Familiarity -0.567 0.003 0.001
## Big5 -0.447 0.003 -0.001 -0.092
## Expressblty -0.387 0.004 0.006 0.014 0.004
## TrialNumber -0.288 -0.027 -0.044 0.006 0.005 -0.012
## Modaltygstr -0.332 -0.010 -0.018 0.004 0.010 0.351 -0.008
## Modalityvcl -0.325 -0.006 -0.014 0.005 -0.003 0.380 -0.047 0.562The Betas seem quite close to how we create the synthetic data, probably there are some other moderations since the causal relations between the variables are quite complex
Open question: dyad + participant
How/whether to include participant in addition to dyad in the random effects
Our thoughts: participant could maybe be nested within dyad, or do we not need participant? Must be something suggested in the literature.
Participant probably must be included because the Big 5 personality trait measure corresponds to each participant within the dyad.
Onur and Frede on nesting:
First, we do not use mono- versus bilingualism or even the group variable as a main independent variable in our model. Instead, it is incorporated in a nested random effect. This indicates that we do not view the speaker group status as a main driver of DM production. Second, the group variable that we utilize does not take one of the three speaker groups as a baseline. We sum-coded the variable which incorporates the idea that monolinguals cannot be an adequate baseline for bilinguals.
Check model diagnostics
performance::check_model(lmer_model)## Failed to compute posterior predictive checks with `re_formula=NULL`.
## Trying again with `re_formula=NA` now.
library(sjPlot)## Warning: Paket 'sjPlot' wurde unter R Version 4.3.3 erstellt## Learn more about sjPlot with 'browseVignettes("sjPlot")'.sjPlot::plot_model(lmer_model)
library(ggplot2)
ggeffects::ggpredict(lmer_model,
terms = c("CommAtt"),
type = "random",
interval = "confidence") |>
as_tibble() |>
ggplot2::ggplot() +
aes(x = x, y = predicted) +
geom_point() +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .2) +
geom_line(group = 1) +
theme_bw()
Going Bayesian
Potential ways:
- LMMs - but we assume there might be some non-linearity (e.g., effort will increase in 2nd attempt, but not third)
- b-splines - but our ‘timeline’ consists of only three points (attempts). We could try but maybe b-splines are unncessarilly complex
- ordinal regression
- see https://solomonkurz.netlify.app/blog/2023-05-21-causal-inference-with-ordinal-regression/ and Chapters 12-13 in McElreath
However, after consultation with Daniela, we can get rid of the non-linear issue by contrast coding
library(ggplot2)
library(patchwork)## Warning: Paket 'patchwork' wurde unter R Version 4.3.3 erstelltlibrary(bayesplot)## Warning: Paket 'bayesplot' wurde unter R Version 4.3.3 erstellt## This is bayesplot version 1.11.1## - Online documentation and vignettes at mc-stan.org/bayesplot## - bayesplot theme set to bayesplot::theme_default()## * Does _not_ affect other ggplot2 plots## * See ?bayesplot_theme_set for details on theme settinglibrary(brms)## Warning: Paket 'brms' wurde unter R Version 4.3.3 erstellt## Lade nötiges Paket: Rcpp## Warning: Paket 'Rcpp' wurde unter R Version 4.3.3 erstellt## Loading 'brms' package (version 2.22.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attache Paket: 'brms'## Das folgende Objekt ist maskiert 'package:bayesplot':
##
## rhat## Das folgende Objekt ist maskiert 'package:lme4':
##
## ngrps## Das folgende Objekt ist maskiert 'package:stats':
##
## arlibrary(beepr)## Warning: Paket 'beepr' wurde unter R Version 4.3.3 erstellt# use all available cores for parallel computing
options(mc.cores = parallel::detectCores())Model 1: starter
Let’s start slowly, first with simple model accounting only for outcome + main predictor
Eff ~ Normal(μ,σ^2) μ= β0 + β1*CommAtt
fit_eff1 <- brm(Eff ~ 1 + CommAtt, data=final_data)
saveRDS(fit_eff1, here("06_Modeling", "models", "fit_eff1.rds"))
beep(5)fit_eff1 <- readRDS(here("06_Modeling", "models", "fit_eff1.rds"))
summary(fit_eff1)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt
## Data: final_data (Number of observations: 5055)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.74 0.01 2.72 2.77 1.00 4017 3144
## CommAtt2M1 1.58 0.03 1.52 1.63 1.00 3578 3162
## CommAtt3M2 -2.12 0.04 -2.19 -2.04 1.00 3399 2958
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.89 0.01 0.88 0.91 1.00 4138 3305
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# even so simple model can quite well capture the relationship in the simulated data
plot(conditional_effects(fit_eff1), points = TRUE)
# but there is quite a lot of variation, causing overlap
pp_check(fit_eff1, type = "dens_overlay")## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# we see that the simulated data are not that loyal to the data
pp_check(fit_eff1, type = "error_scatter_avg")## Using all posterior draws for ppc type 'error_scatter_avg' by default.
# instead of cloud of points, we see almost perfectly linear relationship between outcome and residuals, indicating strong problems with the independence assumption of the errors.Model 2: varying intercept (random effect)
Now we know we will probably need partial pooling for participants, because they share lot of commonalities (Later we might adapt this, because it might be participant within dyads. Dyads might also share lot of commonalities)
fit_eff2 <- brm(Eff ~ 1 + CommAtt + (1 | Participant), data=final_data)
saveRDS(fit_eff2, here("06_Modeling", "models", "fit_eff2.rds"))
beep(5)fit_eff2 <- readRDS(here("06_Modeling", "models", "fit_eff2.rds"))
summary(fit_eff2)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + (1 | Participant)
## Data: final_data (Number of observations: 5055)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.51 0.03 0.45 0.58 1.01 406 795
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.75 0.05 2.65 2.84 1.02 204 383
## CommAtt2M1 1.59 0.02 1.54 1.63 1.00 4096 3460
## CommAtt3M2 -2.12 0.03 -2.18 -2.06 1.00 3988 3397
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.74 0.01 0.72 0.75 1.00 4894 3154
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(conditional_effects(fit_eff2), points = TRUE)
# stays the same
pp_check(fit_eff2, type = "dens_overlay")## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# not really improvement
pp_check(fit_eff2, type = "error_scatter_avg")## Using all posterior draws for ppc type 'error_scatter_avg' by default.
# ok, here we see some improvement here!
# checking to coeff of the parameters
coef(fit_eff2)## $Participant
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 1 3.579610 0.11441565 3.352994 3.812583
## 2 2.522447 0.11044733 2.314125 2.741480
## 3 3.146291 0.10402491 2.944224 3.347758
## 4 2.262994 0.10891029 2.051830 2.477337
## 5 1.668679 0.11156294 1.453740 1.885679
## 6 3.210786 0.10755374 3.001736 3.421982
## 7 1.743726 0.09956315 1.546578 1.942408
## 8 2.415176 0.11570030 2.188098 2.644443
## 9 2.941818 0.10847545 2.726732 3.153139
## 10 2.850522 0.11247866 2.631404 3.071326
## 11 3.235421 0.10197221 3.035279 3.431147
## 12 3.268229 0.10853672 3.055044 3.486600
## 13 3.101527 0.11377116 2.872792 3.323975
## 14 2.134733 0.10677733 1.927514 2.340458
## 15 2.820435 0.11210555 2.603776 3.041106
## 16 2.327921 0.10479141 2.120863 2.528586
## 17 2.612322 0.10537909 2.405861 2.819664
## 18 3.023965 0.11462556 2.795903 3.246746
## 19 2.778095 0.11530711 2.551082 3.002631
## 20 3.151138 0.10798837 2.942186 3.357480
## 21 3.567949 0.11302305 3.352112 3.786573
## 22 2.429940 0.11628572 2.196410 2.659209
## 23 2.968761 0.11250622 2.749033 3.188065
## 24 2.704271 0.10887674 2.491699 2.912234
## 25 2.672293 0.10259883 2.471873 2.872144
## 26 2.779301 0.10857503 2.573211 2.993449
## 27 2.465089 0.11036437 2.249112 2.681764
## 28 2.223275 0.11764174 1.989483 2.456018
## 29 1.846008 0.11587214 1.628159 2.072764
## 30 3.048777 0.11053537 2.831188 3.261401
## 31 2.809188 0.10616442 2.601615 3.014372
## 32 2.543952 0.10470197 2.338840 2.751140
## 33 2.552387 0.11051521 2.331680 2.763790
## 34 2.076140 0.11182684 1.856244 2.292397
## 35 2.219996 0.11277154 1.995945 2.436680
## 36 2.493701 0.10603293 2.281718 2.702046
## 37 3.679949 0.11266136 3.461175 3.900210
## 38 2.736239 0.11295006 2.512656 2.950943
## 39 2.702837 0.11809388 2.469696 2.933019
## 40 2.882159 0.10722781 2.678426 3.093300
## 41 2.472661 0.11187430 2.253757 2.687813
## 42 2.745807 0.11707738 2.510145 2.974322
## 43 2.702553 0.10588311 2.497378 2.912475
## 44 2.260898 0.11366341 2.040076 2.493751
## 45 3.384675 0.11866167 3.152057 3.616854
## 46 2.716896 0.10449897 2.511432 2.925632
## 47 1.852768 0.10364244 1.652043 2.051654
## 48 3.342887 0.11549076 3.118391 3.559941
## 49 1.771033 0.11031665 1.555702 1.992023
## 50 3.346107 0.12714312 3.099476 3.597375
## 51 3.877450 0.11178799 3.661954 4.090479
## 52 1.829994 0.10348166 1.624729 2.030700
## 53 2.396591 0.11912672 2.162229 2.638269
## 54 2.380222 0.11579929 2.150989 2.602008
## 55 3.182091 0.10873405 2.974668 3.395113
## 56 2.315297 0.11311378 2.096352 2.537675
## 57 2.064963 0.11308449 1.837772 2.284046
## 58 2.314193 0.11107981 2.098013 2.523934
## 59 2.399026 0.11424194 2.180342 2.626832
## 60 2.089906 0.11778003 1.866431 2.316164
## 61 2.700226 0.11952011 2.469309 2.926424
## 62 2.619516 0.12051327 2.384422 2.857954
## 63 3.171361 0.11073947 2.951481 3.382997
## 64 2.611329 0.11696504 2.383084 2.841700
## 65 2.688010 0.11691878 2.461584 2.922726
## 66 3.127317 0.10673919 2.919742 3.336950
## 67 2.255830 0.10563610 2.052473 2.461304
## 68 2.857681 0.12025854 2.624437 3.091862
## 69 3.289623 0.11099989 3.076365 3.506857
## 70 2.945625 0.10547291 2.740275 3.152410
## 71 3.512916 0.10935253 3.291953 3.730189
## 72 2.990230 0.10685282 2.780816 3.191803
## 73 1.771854 0.10382751 1.571826 1.971438
## 74 2.756099 0.11608226 2.531082 2.981185
## 75 3.400385 0.11514800 3.176591 3.622344
## 76 2.866797 0.11839811 2.632107 3.091283
## 77 2.640661 0.11011714 2.426721 2.854269
## 78 2.335450 0.10309816 2.133630 2.532908
## 79 3.561743 0.10861026 3.351745 3.772179
## 80 3.247257 0.11588270 3.016624 3.480770
## 81 2.780556 0.10592433 2.576949 2.991000
## 82 1.933096 0.12234709 1.691496 2.181393
## 83 2.621529 0.11831408 2.391282 2.850984
## 84 3.181666 0.11442281 2.960273 3.404832
## 85 3.313460 0.10981852 3.101184 3.523723
## 86 3.107441 0.10816763 2.899903 3.327315
## 87 2.968963 0.10515426 2.768879 3.171367
## 88 2.946239 0.11383147 2.728444 3.167263
## 89 2.801226 0.10178248 2.600620 2.999539
## 90 3.110846 0.11296329 2.897673 3.326217
## 91 2.349701 0.11047740 2.129544 2.561816
## 92 3.065909 0.11395658 2.840459 3.290479
## 93 3.123400 0.11582261 2.903462 3.350258
## 94 2.696007 0.12183158 2.459617 2.932619
## 95 2.789953 0.10861537 2.574905 3.002887
## 96 3.029803 0.10494297 2.822342 3.237188
## 97 3.239874 0.10240338 3.044348 3.440797
## 98 2.175207 0.10657819 1.972430 2.380844
## 99 3.136929 0.11506284 2.914134 3.358832
## 100 1.889527 0.10989286 1.671083 2.102376
## 101 3.043525 0.11398017 2.819380 3.263130
## 102 2.275650 0.10904142 2.071340 2.491571
## 103 2.072891 0.10972387 1.856593 2.286583
## 104 2.762111 0.10780242 2.551544 2.978329
## 105 1.945242 0.11819705 1.708879 2.172253
## 106 3.111959 0.12023253 2.879270 3.347590
## 107 2.747825 0.11557633 2.528067 2.970314
## 108 2.987815 0.11518119 2.764853 3.213802
## 109 3.613064 0.11517567 3.383850 3.836056
## 110 3.726170 0.11354328 3.502244 3.946415
## 111 2.775667 0.11775600 2.542527 3.008275
## 112 2.727688 0.12697156 2.477013 2.969186
## 113 3.086444 0.10924760 2.868714 3.301903
## 114 3.506301 0.11439792 3.281420 3.723125
## 115 2.603470 0.11288770 2.385401 2.817476
## 116 2.761380 0.11791953 2.526261 2.993304
## 117 3.276696 0.10850088 3.067212 3.495088
## 118 2.783191 0.11433916 2.558433 3.005766
## 119 3.048421 0.10547291 2.835668 3.258379
## 120 1.960899 0.10806763 1.747192 2.173225
##
## , , CommAtt2M1
##
## Estimate Est.Error Q2.5 Q97.5
## 1 1.588552 0.02323773 1.542844 1.633812
## 2 1.588552 0.02323773 1.542844 1.633812
## 3 1.588552 0.02323773 1.542844 1.633812
## 4 1.588552 0.02323773 1.542844 1.633812
## 5 1.588552 0.02323773 1.542844 1.633812
## 6 1.588552 0.02323773 1.542844 1.633812
## 7 1.588552 0.02323773 1.542844 1.633812
## 8 1.588552 0.02323773 1.542844 1.633812
## 9 1.588552 0.02323773 1.542844 1.633812
## 10 1.588552 0.02323773 1.542844 1.633812
## 11 1.588552 0.02323773 1.542844 1.633812
## 12 1.588552 0.02323773 1.542844 1.633812
## 13 1.588552 0.02323773 1.542844 1.633812
## 14 1.588552 0.02323773 1.542844 1.633812
## 15 1.588552 0.02323773 1.542844 1.633812
## 16 1.588552 0.02323773 1.542844 1.633812
## 17 1.588552 0.02323773 1.542844 1.633812
## 18 1.588552 0.02323773 1.542844 1.633812
## 19 1.588552 0.02323773 1.542844 1.633812
## 20 1.588552 0.02323773 1.542844 1.633812
## 21 1.588552 0.02323773 1.542844 1.633812
## 22 1.588552 0.02323773 1.542844 1.633812
## 23 1.588552 0.02323773 1.542844 1.633812
## 24 1.588552 0.02323773 1.542844 1.633812
## 25 1.588552 0.02323773 1.542844 1.633812
## 26 1.588552 0.02323773 1.542844 1.633812
## 27 1.588552 0.02323773 1.542844 1.633812
## 28 1.588552 0.02323773 1.542844 1.633812
## 29 1.588552 0.02323773 1.542844 1.633812
## 30 1.588552 0.02323773 1.542844 1.633812
## 31 1.588552 0.02323773 1.542844 1.633812
## 32 1.588552 0.02323773 1.542844 1.633812
## 33 1.588552 0.02323773 1.542844 1.633812
## 34 1.588552 0.02323773 1.542844 1.633812
## 35 1.588552 0.02323773 1.542844 1.633812
## 36 1.588552 0.02323773 1.542844 1.633812
## 37 1.588552 0.02323773 1.542844 1.633812
## 38 1.588552 0.02323773 1.542844 1.633812
## 39 1.588552 0.02323773 1.542844 1.633812
## 40 1.588552 0.02323773 1.542844 1.633812
## 41 1.588552 0.02323773 1.542844 1.633812
## 42 1.588552 0.02323773 1.542844 1.633812
## 43 1.588552 0.02323773 1.542844 1.633812
## 44 1.588552 0.02323773 1.542844 1.633812
## 45 1.588552 0.02323773 1.542844 1.633812
## 46 1.588552 0.02323773 1.542844 1.633812
## 47 1.588552 0.02323773 1.542844 1.633812
## 48 1.588552 0.02323773 1.542844 1.633812
## 49 1.588552 0.02323773 1.542844 1.633812
## 50 1.588552 0.02323773 1.542844 1.633812
## 51 1.588552 0.02323773 1.542844 1.633812
## 52 1.588552 0.02323773 1.542844 1.633812
## 53 1.588552 0.02323773 1.542844 1.633812
## 54 1.588552 0.02323773 1.542844 1.633812
## 55 1.588552 0.02323773 1.542844 1.633812
## 56 1.588552 0.02323773 1.542844 1.633812
## 57 1.588552 0.02323773 1.542844 1.633812
## 58 1.588552 0.02323773 1.542844 1.633812
## 59 1.588552 0.02323773 1.542844 1.633812
## 60 1.588552 0.02323773 1.542844 1.633812
## 61 1.588552 0.02323773 1.542844 1.633812
## 62 1.588552 0.02323773 1.542844 1.633812
## 63 1.588552 0.02323773 1.542844 1.633812
## 64 1.588552 0.02323773 1.542844 1.633812
## 65 1.588552 0.02323773 1.542844 1.633812
## 66 1.588552 0.02323773 1.542844 1.633812
## 67 1.588552 0.02323773 1.542844 1.633812
## 68 1.588552 0.02323773 1.542844 1.633812
## 69 1.588552 0.02323773 1.542844 1.633812
## 70 1.588552 0.02323773 1.542844 1.633812
## 71 1.588552 0.02323773 1.542844 1.633812
## 72 1.588552 0.02323773 1.542844 1.633812
## 73 1.588552 0.02323773 1.542844 1.633812
## 74 1.588552 0.02323773 1.542844 1.633812
## 75 1.588552 0.02323773 1.542844 1.633812
## 76 1.588552 0.02323773 1.542844 1.633812
## 77 1.588552 0.02323773 1.542844 1.633812
## 78 1.588552 0.02323773 1.542844 1.633812
## 79 1.588552 0.02323773 1.542844 1.633812
## 80 1.588552 0.02323773 1.542844 1.633812
## 81 1.588552 0.02323773 1.542844 1.633812
## 82 1.588552 0.02323773 1.542844 1.633812
## 83 1.588552 0.02323773 1.542844 1.633812
## 84 1.588552 0.02323773 1.542844 1.633812
## 85 1.588552 0.02323773 1.542844 1.633812
## 86 1.588552 0.02323773 1.542844 1.633812
## 87 1.588552 0.02323773 1.542844 1.633812
## 88 1.588552 0.02323773 1.542844 1.633812
## 89 1.588552 0.02323773 1.542844 1.633812
## 90 1.588552 0.02323773 1.542844 1.633812
## 91 1.588552 0.02323773 1.542844 1.633812
## 92 1.588552 0.02323773 1.542844 1.633812
## 93 1.588552 0.02323773 1.542844 1.633812
## 94 1.588552 0.02323773 1.542844 1.633812
## 95 1.588552 0.02323773 1.542844 1.633812
## 96 1.588552 0.02323773 1.542844 1.633812
## 97 1.588552 0.02323773 1.542844 1.633812
## 98 1.588552 0.02323773 1.542844 1.633812
## 99 1.588552 0.02323773 1.542844 1.633812
## 100 1.588552 0.02323773 1.542844 1.633812
## 101 1.588552 0.02323773 1.542844 1.633812
## 102 1.588552 0.02323773 1.542844 1.633812
## 103 1.588552 0.02323773 1.542844 1.633812
## 104 1.588552 0.02323773 1.542844 1.633812
## 105 1.588552 0.02323773 1.542844 1.633812
## 106 1.588552 0.02323773 1.542844 1.633812
## 107 1.588552 0.02323773 1.542844 1.633812
## 108 1.588552 0.02323773 1.542844 1.633812
## 109 1.588552 0.02323773 1.542844 1.633812
## 110 1.588552 0.02323773 1.542844 1.633812
## 111 1.588552 0.02323773 1.542844 1.633812
## 112 1.588552 0.02323773 1.542844 1.633812
## 113 1.588552 0.02323773 1.542844 1.633812
## 114 1.588552 0.02323773 1.542844 1.633812
## 115 1.588552 0.02323773 1.542844 1.633812
## 116 1.588552 0.02323773 1.542844 1.633812
## 117 1.588552 0.02323773 1.542844 1.633812
## 118 1.588552 0.02323773 1.542844 1.633812
## 119 1.588552 0.02323773 1.542844 1.633812
## 120 1.588552 0.02323773 1.542844 1.633812
##
## , , CommAtt3M2
##
## Estimate Est.Error Q2.5 Q97.5
## 1 -2.120611 0.03038427 -2.180229 -2.061068
## 2 -2.120611 0.03038427 -2.180229 -2.061068
## 3 -2.120611 0.03038427 -2.180229 -2.061068
## 4 -2.120611 0.03038427 -2.180229 -2.061068
## 5 -2.120611 0.03038427 -2.180229 -2.061068
## 6 -2.120611 0.03038427 -2.180229 -2.061068
## 7 -2.120611 0.03038427 -2.180229 -2.061068
## 8 -2.120611 0.03038427 -2.180229 -2.061068
## 9 -2.120611 0.03038427 -2.180229 -2.061068
## 10 -2.120611 0.03038427 -2.180229 -2.061068
## 11 -2.120611 0.03038427 -2.180229 -2.061068
## 12 -2.120611 0.03038427 -2.180229 -2.061068
## 13 -2.120611 0.03038427 -2.180229 -2.061068
## 14 -2.120611 0.03038427 -2.180229 -2.061068
## 15 -2.120611 0.03038427 -2.180229 -2.061068
## 16 -2.120611 0.03038427 -2.180229 -2.061068
## 17 -2.120611 0.03038427 -2.180229 -2.061068
## 18 -2.120611 0.03038427 -2.180229 -2.061068
## 19 -2.120611 0.03038427 -2.180229 -2.061068
## 20 -2.120611 0.03038427 -2.180229 -2.061068
## 21 -2.120611 0.03038427 -2.180229 -2.061068
## 22 -2.120611 0.03038427 -2.180229 -2.061068
## 23 -2.120611 0.03038427 -2.180229 -2.061068
## 24 -2.120611 0.03038427 -2.180229 -2.061068
## 25 -2.120611 0.03038427 -2.180229 -2.061068
## 26 -2.120611 0.03038427 -2.180229 -2.061068
## 27 -2.120611 0.03038427 -2.180229 -2.061068
## 28 -2.120611 0.03038427 -2.180229 -2.061068
## 29 -2.120611 0.03038427 -2.180229 -2.061068
## 30 -2.120611 0.03038427 -2.180229 -2.061068
## 31 -2.120611 0.03038427 -2.180229 -2.061068
## 32 -2.120611 0.03038427 -2.180229 -2.061068
## 33 -2.120611 0.03038427 -2.180229 -2.061068
## 34 -2.120611 0.03038427 -2.180229 -2.061068
## 35 -2.120611 0.03038427 -2.180229 -2.061068
## 36 -2.120611 0.03038427 -2.180229 -2.061068
## 37 -2.120611 0.03038427 -2.180229 -2.061068
## 38 -2.120611 0.03038427 -2.180229 -2.061068
## 39 -2.120611 0.03038427 -2.180229 -2.061068
## 40 -2.120611 0.03038427 -2.180229 -2.061068
## 41 -2.120611 0.03038427 -2.180229 -2.061068
## 42 -2.120611 0.03038427 -2.180229 -2.061068
## 43 -2.120611 0.03038427 -2.180229 -2.061068
## 44 -2.120611 0.03038427 -2.180229 -2.061068
## 45 -2.120611 0.03038427 -2.180229 -2.061068
## 46 -2.120611 0.03038427 -2.180229 -2.061068
## 47 -2.120611 0.03038427 -2.180229 -2.061068
## 48 -2.120611 0.03038427 -2.180229 -2.061068
## 49 -2.120611 0.03038427 -2.180229 -2.061068
## 50 -2.120611 0.03038427 -2.180229 -2.061068
## 51 -2.120611 0.03038427 -2.180229 -2.061068
## 52 -2.120611 0.03038427 -2.180229 -2.061068
## 53 -2.120611 0.03038427 -2.180229 -2.061068
## 54 -2.120611 0.03038427 -2.180229 -2.061068
## 55 -2.120611 0.03038427 -2.180229 -2.061068
## 56 -2.120611 0.03038427 -2.180229 -2.061068
## 57 -2.120611 0.03038427 -2.180229 -2.061068
## 58 -2.120611 0.03038427 -2.180229 -2.061068
## 59 -2.120611 0.03038427 -2.180229 -2.061068
## 60 -2.120611 0.03038427 -2.180229 -2.061068
## 61 -2.120611 0.03038427 -2.180229 -2.061068
## 62 -2.120611 0.03038427 -2.180229 -2.061068
## 63 -2.120611 0.03038427 -2.180229 -2.061068
## 64 -2.120611 0.03038427 -2.180229 -2.061068
## 65 -2.120611 0.03038427 -2.180229 -2.061068
## 66 -2.120611 0.03038427 -2.180229 -2.061068
## 67 -2.120611 0.03038427 -2.180229 -2.061068
## 68 -2.120611 0.03038427 -2.180229 -2.061068
## 69 -2.120611 0.03038427 -2.180229 -2.061068
## 70 -2.120611 0.03038427 -2.180229 -2.061068
## 71 -2.120611 0.03038427 -2.180229 -2.061068
## 72 -2.120611 0.03038427 -2.180229 -2.061068
## 73 -2.120611 0.03038427 -2.180229 -2.061068
## 74 -2.120611 0.03038427 -2.180229 -2.061068
## 75 -2.120611 0.03038427 -2.180229 -2.061068
## 76 -2.120611 0.03038427 -2.180229 -2.061068
## 77 -2.120611 0.03038427 -2.180229 -2.061068
## 78 -2.120611 0.03038427 -2.180229 -2.061068
## 79 -2.120611 0.03038427 -2.180229 -2.061068
## 80 -2.120611 0.03038427 -2.180229 -2.061068
## 81 -2.120611 0.03038427 -2.180229 -2.061068
## 82 -2.120611 0.03038427 -2.180229 -2.061068
## 83 -2.120611 0.03038427 -2.180229 -2.061068
## 84 -2.120611 0.03038427 -2.180229 -2.061068
## 85 -2.120611 0.03038427 -2.180229 -2.061068
## 86 -2.120611 0.03038427 -2.180229 -2.061068
## 87 -2.120611 0.03038427 -2.180229 -2.061068
## 88 -2.120611 0.03038427 -2.180229 -2.061068
## 89 -2.120611 0.03038427 -2.180229 -2.061068
## 90 -2.120611 0.03038427 -2.180229 -2.061068
## 91 -2.120611 0.03038427 -2.180229 -2.061068
## 92 -2.120611 0.03038427 -2.180229 -2.061068
## 93 -2.120611 0.03038427 -2.180229 -2.061068
## 94 -2.120611 0.03038427 -2.180229 -2.061068
## 95 -2.120611 0.03038427 -2.180229 -2.061068
## 96 -2.120611 0.03038427 -2.180229 -2.061068
## 97 -2.120611 0.03038427 -2.180229 -2.061068
## 98 -2.120611 0.03038427 -2.180229 -2.061068
## 99 -2.120611 0.03038427 -2.180229 -2.061068
## 100 -2.120611 0.03038427 -2.180229 -2.061068
## 101 -2.120611 0.03038427 -2.180229 -2.061068
## 102 -2.120611 0.03038427 -2.180229 -2.061068
## 103 -2.120611 0.03038427 -2.180229 -2.061068
## 104 -2.120611 0.03038427 -2.180229 -2.061068
## 105 -2.120611 0.03038427 -2.180229 -2.061068
## 106 -2.120611 0.03038427 -2.180229 -2.061068
## 107 -2.120611 0.03038427 -2.180229 -2.061068
## 108 -2.120611 0.03038427 -2.180229 -2.061068
## 109 -2.120611 0.03038427 -2.180229 -2.061068
## 110 -2.120611 0.03038427 -2.180229 -2.061068
## 111 -2.120611 0.03038427 -2.180229 -2.061068
## 112 -2.120611 0.03038427 -2.180229 -2.061068
## 113 -2.120611 0.03038427 -2.180229 -2.061068
## 114 -2.120611 0.03038427 -2.180229 -2.061068
## 115 -2.120611 0.03038427 -2.180229 -2.061068
## 116 -2.120611 0.03038427 -2.180229 -2.061068
## 117 -2.120611 0.03038427 -2.180229 -2.061068
## 118 -2.120611 0.03038427 -2.180229 -2.061068
## 119 -2.120611 0.03038427 -2.180229 -2.061068
## 120 -2.120611 0.03038427 -2.180229 -2.061068Some more diagnostics
vars <- c("b_Intercept", "r_Participant[1,Intercept]", "r_Participant[2,Intercept]")
pairs(fit_eff2, variable = vars)
We see negative correlation between b0 and b0j (and b0j have weak
positive correlation with one another)
draws_eff2 <- as.matrix(fit_eff2, variable = vars)
colnames(draws_eff2) <- c("b0", "b0_1", "b0_2")
cors_eff2 <- cor(draws_eff2)
print(cors_eff2, digits = 2)## b0 b0_1 b0_2
## b0 1.00 -0.36 -0.39
## b0_1 -0.36 1.00 0.12
## b0_2 -0.39 0.12 1.00Correlation between the parameters
coef_eff2 <- coef(fit_eff2, summary = FALSE)
print(cor(coef_eff2$Participant[, 1:5, "Intercept"]), digits = 2)## 1 2 3 4 5
## 1 1.000 -0.0178 0.0461 0.012 -0.0100
## 2 -0.018 1.0000 0.0039 0.063 -0.0084
## 3 0.046 0.0039 1.0000 0.057 0.0135
## 4 0.012 0.0626 0.0571 1.000 0.0378
## 5 -0.010 -0.0084 0.0135 0.038 1.0000Per-subject predicitons into the plot
conditions <- make_conditions(final_data, "Participant")
ce <- conditional_effects(
fit_eff2, conditions = conditions,
re_formula = NULL
)
plot(ce, ncol = 6, points = TRUE)
# ok this is too many pcns to see anything:)Compare the first two models
loo_eff1 <- loo(fit_eff1)
loo_eff2 <- loo(fit_eff2)
loo_compare(loo_eff1, loo_eff2)## elpd_diff se_diff
## fit_eff2 0.0 0.0
## fit_eff1 -912.1 38.1Model 3: varying slope
Let’s also add varying slope, because ComAtt will too vary across subjects
fit_eff3 <- brm(Eff ~ 1 + CommAtt + (1 + CommAtt | Participant), # note that there was a mistake before
data=final_data,
# Number of cores to use
cores = 4,)
saveRDS(fit_eff3, here("06_Modeling", "models", "fit_eff3.rds"))
beep(5)fit_eff3 <- readRDS(here("06_Modeling", "models", "fit_eff3.rds"))
summary(fit_eff3)## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.## Warning: There were 1278 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + (1 | CommAtt + Participant)
## Data: final_data (Number of observations: 5055)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~CommAtt (Number of levels: 3)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.72 1.94 0.08 5.71 1.27 11 55
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.53 0.03 0.45 0.59 1.16 19 436
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.23 1.69 0.05 6.82 1.50 7 32
## CommAtt2M1 2.10 2.46 -3.39 6.78 1.15 61 668
## CommAtt3M2 -1.25 2.86 -7.81 3.37 1.14 20 106
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.74 0.01 0.72 0.75 1.41 557 1789
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(conditional_effects(fit_eff3), points = TRUE)
# now the condidence intervals are much wider
pp_check(fit_eff3, type = "dens_overlay")## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# not really improvement
pp_check(fit_eff3, type = "error_scatter_avg")## Using all posterior draws for ppc type 'error_scatter_avg' by default.
# blobby, but still very correlatedModel 4: full DAG
Now let’s finally replicate our DAG
fit_eff4 <- brm(Eff ~ 1 + CommAtt + Familiarity + Big5 + Expressibility_z + TrialNumber_c + Modality + (1 | Participant) + (1 | Concept),
data = final_data,
cores = 4)
saveRDS(fit_eff4, here("06_Modeling", "models", "fit_eff4.rds"))
beep(5)fit_eff4 <- readRDS(here("06_Modeling", "models", "fit_eff4.rds"))
# summary
summary(fit_eff4)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Familiarity + Big5 + Expressibility_z + TrialNumber_c + Modality + (1 | Participant) + (1 | Concept)
## Data: final_data (Number of observations: 5055)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.03 0.01 0.00 0.06 1.00 1309 1353
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.02 0.01 0.00 0.04 1.00 1886 2247
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.71 0.01 2.69 2.74 1.00 4845 3442
## CommAtt2M1 1.60 0.02 1.56 1.64 1.00 6017 3209
## CommAtt3M2 -2.12 0.03 -2.17 -2.07 1.00 6121 3464
## Familiarity 1.23 0.03 1.16 1.29 1.00 5385 3222
## Big5 1.19 0.03 1.13 1.25 1.00 5173 3189
## Expressibility_z 0.39 0.01 0.37 0.41 1.00 4850 2713
## TrialNumber_c -0.00 0.00 -0.00 0.00 1.00 3868 2355
## Modality1 -0.55 0.03 -0.60 -0.49 1.00 5645 2824
## Modality2 0.30 0.03 0.25 0.35 1.00 5678 3188
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.65 0.01 0.63 0.66 1.00 6786 2916
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# the coefficients again look similar to what we set in the data (except perhaps expressibility?)
plot(fit_eff4)
plot(conditional_effects(fit_eff4), points = TRUE)
# the CrI are now again in sensible width
pp_check(fit_eff4, type = "dens_overlay")## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# not really improvement
pp_check(fit_eff4, type = "error_scatter_avg")## Using all posterior draws for ppc type 'error_scatter_avg' by default.
# blobby, but still correlated
# note that we are still missing proper non-def priorsModel 5: full + varying intercepts
fit_eff5 <- brm(Eff ~ 1 + CommAtt + Modality + (1 | Participant) + (1 | Concept) + (1 | Familiarity) + (1 | Big5) + (1 | Expressibility_z) + (1 | TrialNumber_c),
data = final_data,
cores = 4)
# NOTE: I accidentally rewrote fit5 with model fit6 so this will need to be re-run and re-saved
saveRDS(fit_eff5, here("06_Modeling", "models", "fit_eff5.rds"))
beep(5)fit_eff5 <- readRDS(here("06_Modeling", "models", "fit_eff5.rds"))
# summary
summary(fit_eff5)## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.## Warning: There were 41 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c)
## Data: final_data (Number of observations: 5055)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Big5 (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.06 0.05 0.00 0.17 1.02 140
## sd(CommAtt2M1) 0.03 0.03 0.00 0.10 1.02 403
## sd(CommAtt3M2) 0.03 0.03 0.00 0.11 1.01 361
## cor(Intercept,CommAtt2M1) 0.03 0.51 -0.87 0.91 1.01 874
## cor(Intercept,CommAtt3M2) -0.15 0.52 -0.94 0.83 1.01 616
## cor(CommAtt2M1,CommAtt3M2) -0.18 0.52 -0.95 0.84 1.01 632
## Tail_ESS
## sd(Intercept) 104
## sd(CommAtt2M1) 151
## sd(CommAtt3M2) 131
## cor(Intercept,CommAtt2M1) 1145
## cor(Intercept,CommAtt3M2) 443
## cor(CommAtt2M1,CommAtt3M2) 286
##
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.06 0.04 0.00 0.13 1.01 338 1047
##
## ~Expressibility_z (Number of levels: 63)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.49 0.07 0.38 0.63 1.03 170
## sd(CommAtt2M1) 0.61 0.06 0.51 0.72 1.01 624
## sd(CommAtt3M2) 0.35 0.04 0.28 0.43 1.01 779
## cor(Intercept,CommAtt2M1) 0.92 0.06 0.79 0.98 1.02 128
## cor(Intercept,CommAtt3M2) -0.95 0.04 -1.00 -0.86 1.01 294
## cor(CommAtt2M1,CommAtt3M2) -0.81 0.06 -0.91 -0.68 1.00 1107
## Tail_ESS
## sd(Intercept) 261
## sd(CommAtt2M1) 1230
## sd(CommAtt3M2) 1843
## cor(Intercept,CommAtt2M1) 172
## cor(Intercept,CommAtt3M2) 614
## cor(CommAtt2M1,CommAtt3M2) 1761
##
## ~Familiarity (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.16 0.18 0.00 0.51 1.56 7
## sd(CommAtt2M1) 0.09 0.11 0.00 0.32 1.53 7
## sd(CommAtt3M2) 0.12 0.15 0.00 0.43 1.53 7
## cor(Intercept,CommAtt2M1) 0.26 0.59 -0.83 0.99 1.50 7
## cor(Intercept,CommAtt3M2) -0.38 0.55 -1.00 0.78 1.51 7
## cor(CommAtt2M1,CommAtt3M2) -0.40 0.55 -1.00 0.77 1.51 7
## Tail_ESS
## sd(Intercept) 39
## sd(CommAtt2M1) 29
## sd(CommAtt3M2) 29
## cor(Intercept,CommAtt2M1) 39
## cor(Intercept,CommAtt3M2) 39
## cor(CommAtt2M1,CommAtt3M2) 32
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.36 0.18 0.01 0.53 1.54 7
## sd(CommAtt2M1) 0.23 0.12 0.00 0.34 1.53 7
## sd(CommAtt3M2) 0.30 0.16 0.00 0.45 1.53 7
## cor(Intercept,CommAtt2M1) 0.71 0.49 -0.70 0.99 1.51 7
## cor(Intercept,CommAtt3M2) -0.76 0.44 -1.00 0.57 1.52 7
## cor(CommAtt2M1,CommAtt3M2) -0.78 0.43 -1.00 0.57 1.52 7
## Tail_ESS
## sd(Intercept) 27
## sd(CommAtt2M1) 31
## sd(CommAtt3M2) 29
## cor(Intercept,CommAtt2M1) 29
## cor(Intercept,CommAtt3M2) 28
## cor(CommAtt2M1,CommAtt3M2) 28
##
## ~TrialNumber_c (Number of levels: 54)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.02 0.01 0.00 0.05 1.00 1085 1283
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.76 0.07 2.61 2.90 1.00 361 650
## CommAtt2M1 1.62 0.08 1.46 1.78 1.00 414 698
## CommAtt3M2 -2.15 0.06 -2.26 -2.02 1.01 417 808
## Modality1 -0.87 0.20 -1.24 -0.49 1.02 129 123
## Modality2 0.45 0.11 0.23 0.66 1.02 190 235
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.56 0.01 0.54 0.57 1.00 3809 1448
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# the coefficients again look similar to what we set in the data (except perhaps expressibility?)
plot(fit_eff5)
plot(conditional_effects(fit_eff5), points = TRUE)
# the CrI are now again in sensible width
pp_check(fit_eff5, type = "dens_overlay")## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# not really improvement
pp_check(fit_eff5, type = "error_scatter_avg")## Using all posterior draws for ppc type 'error_scatter_avg' by default.
# blobby, but still correlatedModel 6: full with varying slopes and intercepts
So we not only assume that different participants have different base level for effort. But also that they differently ‘response’ to the communicative attempt
Same for familiarity, Big5, expressibility
fit_eff6 <- brm(Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c),
data = final_data,
cores = 4)
saveRDS(fit_eff6, here("06_Modeling", "models", "fit_eff6.rds"))
beep(5)fit_eff6 <- readRDS(here("06_Modeling", "models", "fit_eff6.rds"))
# summary
summary(fit_eff6)## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.## Warning: There were 41 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c)
## Data: final_data (Number of observations: 5055)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Big5 (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.06 0.05 0.00 0.17 1.02 140
## sd(CommAtt2M1) 0.03 0.03 0.00 0.10 1.02 403
## sd(CommAtt3M2) 0.03 0.03 0.00 0.11 1.01 361
## cor(Intercept,CommAtt2M1) 0.03 0.51 -0.87 0.91 1.01 874
## cor(Intercept,CommAtt3M2) -0.15 0.52 -0.94 0.83 1.01 616
## cor(CommAtt2M1,CommAtt3M2) -0.18 0.52 -0.95 0.84 1.01 632
## Tail_ESS
## sd(Intercept) 104
## sd(CommAtt2M1) 151
## sd(CommAtt3M2) 131
## cor(Intercept,CommAtt2M1) 1145
## cor(Intercept,CommAtt3M2) 443
## cor(CommAtt2M1,CommAtt3M2) 286
##
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.06 0.04 0.00 0.13 1.01 338 1047
##
## ~Expressibility_z (Number of levels: 63)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.49 0.07 0.38 0.63 1.03 170
## sd(CommAtt2M1) 0.61 0.06 0.51 0.72 1.01 624
## sd(CommAtt3M2) 0.35 0.04 0.28 0.43 1.01 779
## cor(Intercept,CommAtt2M1) 0.92 0.06 0.79 0.98 1.02 128
## cor(Intercept,CommAtt3M2) -0.95 0.04 -1.00 -0.86 1.01 294
## cor(CommAtt2M1,CommAtt3M2) -0.81 0.06 -0.91 -0.68 1.00 1107
## Tail_ESS
## sd(Intercept) 261
## sd(CommAtt2M1) 1230
## sd(CommAtt3M2) 1843
## cor(Intercept,CommAtt2M1) 172
## cor(Intercept,CommAtt3M2) 614
## cor(CommAtt2M1,CommAtt3M2) 1761
##
## ~Familiarity (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.16 0.18 0.00 0.51 1.56 7
## sd(CommAtt2M1) 0.09 0.11 0.00 0.32 1.53 7
## sd(CommAtt3M2) 0.12 0.15 0.00 0.43 1.53 7
## cor(Intercept,CommAtt2M1) 0.26 0.59 -0.83 0.99 1.50 7
## cor(Intercept,CommAtt3M2) -0.38 0.55 -1.00 0.78 1.51 7
## cor(CommAtt2M1,CommAtt3M2) -0.40 0.55 -1.00 0.77 1.51 7
## Tail_ESS
## sd(Intercept) 39
## sd(CommAtt2M1) 29
## sd(CommAtt3M2) 29
## cor(Intercept,CommAtt2M1) 39
## cor(Intercept,CommAtt3M2) 39
## cor(CommAtt2M1,CommAtt3M2) 32
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.36 0.18 0.01 0.53 1.54 7
## sd(CommAtt2M1) 0.23 0.12 0.00 0.34 1.53 7
## sd(CommAtt3M2) 0.30 0.16 0.00 0.45 1.53 7
## cor(Intercept,CommAtt2M1) 0.71 0.49 -0.70 0.99 1.51 7
## cor(Intercept,CommAtt3M2) -0.76 0.44 -1.00 0.57 1.52 7
## cor(CommAtt2M1,CommAtt3M2) -0.78 0.43 -1.00 0.57 1.52 7
## Tail_ESS
## sd(Intercept) 27
## sd(CommAtt2M1) 31
## sd(CommAtt3M2) 29
## cor(Intercept,CommAtt2M1) 29
## cor(Intercept,CommAtt3M2) 28
## cor(CommAtt2M1,CommAtt3M2) 28
##
## ~TrialNumber_c (Number of levels: 54)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.02 0.01 0.00 0.05 1.00 1085 1283
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.76 0.07 2.61 2.90 1.00 361 650
## CommAtt2M1 1.62 0.08 1.46 1.78 1.00 414 698
## CommAtt3M2 -2.15 0.06 -2.26 -2.02 1.01 417 808
## Modality1 -0.87 0.20 -1.24 -0.49 1.02 129 123
## Modality2 0.45 0.11 0.23 0.66 1.02 190 235
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.56 0.01 0.54 0.57 1.00 3809 1448
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# not that there are troubles with convergence! (default priors)
plot(fit_eff6)
plot(conditional_effects(fit_eff6), points = TRUE)
# looks ok
pp_check(fit_eff6, type = "dens_overlay")## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# nice!!! this is definitely improvement
pp_check(fit_eff6, type = "error_scatter_avg")## Using all posterior draws for ppc type 'error_scatter_avg' by default.
# even this looks quite niceModel 7: adding our own priors
Now we will try to fit the same model but with our priors
# check what priors need to be specified
get_prior(Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c),
prior = priors_eff7,
data = final_data)## prior class coef group resp dpar nlpar
## (flat) b
## (flat) b CommAtt2M1
## (flat) b CommAtt3M2
## (flat) b Modality1
## (flat) b Modality2
## lkj(1) cor
## lkj(1) cor Big5
## lkj(1) cor Expressibility_z
## lkj(1) cor Familiarity
## lkj(1) cor Participant
## student_t(3, 2.6, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd Big5
## student_t(3, 0, 2.5) sd CommAtt2M1 Big5
## student_t(3, 0, 2.5) sd CommAtt3M2 Big5
## student_t(3, 0, 2.5) sd Intercept Big5
## student_t(3, 0, 2.5) sd Concept
## student_t(3, 0, 2.5) sd Intercept Concept
## student_t(3, 0, 2.5) sd Expressibility_z
## student_t(3, 0, 2.5) sd CommAtt2M1 Expressibility_z
## student_t(3, 0, 2.5) sd CommAtt3M2 Expressibility_z
## student_t(3, 0, 2.5) sd Intercept Expressibility_z
## student_t(3, 0, 2.5) sd Familiarity
## student_t(3, 0, 2.5) sd CommAtt2M1 Familiarity
## student_t(3, 0, 2.5) sd CommAtt3M2 Familiarity
## student_t(3, 0, 2.5) sd Intercept Familiarity
## student_t(3, 0, 2.5) sd Participant
## student_t(3, 0, 2.5) sd CommAtt2M1 Participant
## student_t(3, 0, 2.5) sd CommAtt3M2 Participant
## student_t(3, 0, 2.5) sd Intercept Participant
## student_t(3, 0, 2.5) sd TrialNumber_c
## student_t(3, 0, 2.5) sd Intercept TrialNumber_c
## student_t(3, 0, 2.5) sigma
## lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 default# Define the priors
priors_eff7 <- c(
prior(normal(2.5, 1), class = "Intercept"), # Prior for the intercept (baseline Effort)
prior(normal(0, 1.5), class = "b", coef = "CommAtt2M1"), # Prior for CommAtt level 2
prior(normal(0, 1.5), class = "b", coef = "CommAtt3M2"), # Prior for CommAtt level 3
prior(normal(0, 0.5), class = "b", coef = "Modality1"), # Prior for Modality level 1
prior(normal(0, 0.5), class = "b", coef = "Modality2"), # Prior for Modality level 2
prior(normal(0, 0.5), class = "sd", group = "Participant"), # Random effect for Participant (Intercept)
prior(normal(0, 0.5), class = "sd", group = "Concept"), # Random effect for Concept (Intercept)
prior(normal(0, 0.5), class = "sd", group = "Familiarity"), # Random effect for Familiarity (Intercept)
prior(normal(0, 0.5), class = "sd", group = "Big5"), # Random effect for Big5 (Intercept)
prior(normal(0, 0.5), class = "sd", group = "Expressibility_z"), # Random effect for Expressibility (Intercept)
prior(normal(0, 0.5), class = "sd", group = "TrialNumber_c"), # Random effect for Trial Number (Intercept)
prior(exponential(1), class = "sigma"), # Prior for residuals
prior(lkj(4), class = "cor")
)# fit model
fit_eff7 <- brm(Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c),
data = final_data,
prior = priors_eff7,
cores = 4)
saveRDS(fit_eff7, here("06_Modeling", "models", "fit_eff7.rds"))
beep(5)fit_eff7 <- readRDS(here("06_Modeling", "models", "fit_eff7.rds"))
# summary
summary(fit_eff7)## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.## Warning: There were 1 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 + CommAtt | Big5) + (1 + CommAtt | Expressibility_z) + (1 | TrialNumber_c)
## Data: final_data (Number of observations: 5055)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Big5 (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.15 0.18 0.00 0.49 1.55 7
## sd(CommAtt2M1) 0.09 0.11 0.00 0.31 1.53 7
## sd(CommAtt3M2) 0.11 0.15 0.00 0.40 1.53 7
## cor(Intercept,CommAtt2M1) 0.23 0.48 -0.59 0.95 1.53 7
## cor(Intercept,CommAtt3M2) -0.27 0.48 -0.97 0.53 1.53 7
## cor(CommAtt2M1,CommAtt3M2) -0.28 0.48 -0.98 0.52 1.53 7
## Tail_ESS
## sd(Intercept) 29
## sd(CommAtt2M1) 29
## sd(CommAtt3M2) 28
## cor(Intercept,CommAtt2M1) 31
## cor(Intercept,CommAtt3M2) 36
## cor(CommAtt2M1,CommAtt3M2) 29
##
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.05 0.03 0.00 0.12 1.01 549 1664
##
## ~Expressibility_z (Number of levels: 63)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.41 0.04 0.33 0.50 1.01 488
## sd(CommAtt2M1) 0.57 0.05 0.48 0.67 1.01 1153
## sd(CommAtt3M2) 0.33 0.03 0.26 0.40 1.01 1652
## cor(Intercept,CommAtt2M1) 0.81 0.09 0.58 0.94 1.01 313
## cor(Intercept,CommAtt3M2) -0.92 0.04 -0.98 -0.82 1.00 1036
## cor(CommAtt2M1,CommAtt3M2) -0.73 0.07 -0.85 -0.57 1.00 1690
## Tail_ESS
## sd(Intercept) 1110
## sd(CommAtt2M1) 1819
## sd(CommAtt3M2) 2925
## cor(Intercept,CommAtt2M1) 320
## cor(Intercept,CommAtt3M2) 2005
## cor(CommAtt2M1,CommAtt3M2) 2606
##
## ~Familiarity (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.15 0.18 0.00 0.49 1.54 7
## sd(CommAtt2M1) 0.09 0.11 0.00 0.31 1.53 7
## sd(CommAtt3M2) 0.11 0.15 0.00 0.41 1.53 7
## cor(Intercept,CommAtt2M1) 0.23 0.47 -0.57 0.95 1.53 7
## cor(Intercept,CommAtt3M2) -0.27 0.48 -0.97 0.54 1.53 7
## cor(CommAtt2M1,CommAtt3M2) -0.28 0.48 -0.98 0.53 1.53 7
## Tail_ESS
## sd(Intercept) 37
## sd(CommAtt2M1) 30
## sd(CommAtt3M2) 30
## cor(Intercept,CommAtt2M1) 31
## cor(Intercept,CommAtt3M2) 29
## cor(CommAtt2M1,CommAtt3M2) 30
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.25 0.21 0.00 0.50 1.74 6
## sd(CommAtt2M1) 0.15 0.13 0.00 0.32 1.73 6
## sd(CommAtt3M2) 0.20 0.17 0.00 0.42 1.73 6
## cor(Intercept,CommAtt2M1) 0.45 0.50 -0.50 0.95 1.73 6
## cor(Intercept,CommAtt3M2) -0.49 0.50 -0.98 0.49 1.73 6
## cor(CommAtt2M1,CommAtt3M2) -0.51 0.50 -0.98 0.47 1.73 6
## Tail_ESS
## sd(Intercept) 168
## sd(CommAtt2M1) 171
## sd(CommAtt3M2) 170
## cor(Intercept,CommAtt2M1) 113
## cor(Intercept,CommAtt3M2) 133
## cor(CommAtt2M1,CommAtt3M2) 125
##
## ~TrialNumber_c (Number of levels: 54)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.02 0.01 0.00 0.04 1.00 1390 1638
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.75 0.07 2.62 2.88 1.02 535 960
## CommAtt2M1 1.61 0.08 1.45 1.75 1.01 666 992
## CommAtt3M2 -2.13 0.06 -2.25 -2.03 1.01 664 1105
## Modality1 -0.64 0.19 -0.98 -0.24 1.01 298 281
## Modality2 0.33 0.11 0.10 0.54 1.01 385 510
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.56 0.01 0.54 0.57 1.00 7878 3132
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# not that there are troubles with convergence! (default priors)
plot(fit_eff7)
# looks scheisse, divergent chains
plot(conditional_effects(fit_eff7), points = TRUE)
# nicht so gut
pp_check(fit_eff7, type = "dens_overlay")## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# stays the same
pp_check(fit_eff7, type = "error_scatter_avg")## Using all posterior draws for ppc type 'error_scatter_avg' by default.
# Note that adding varying slope is maybe not necessary for all the current variables, namely Big5, Expressibility, TrialNumber and Familiarity.
We need to have them in the model to avoid confounds (see DAG), but let’s take them one by one
- Familiarity
As it’s influencing CommAtt and Effort, it should be there as fixed effect to block the causal path.
We also expect different baselines of effort depending on how familiar two people are, so varying intercept is accurate for this
Do we expect that familiarity affects people differently? Perhaps not necessarily. However, for some people it might be that familiarity leads to bigger willingness to be more effortful, for some it means that they know each other so well that they don’t need to put so much effort into things (they are more conventional in a sense)
There should also potentially be an interaction between Familiarity and CommAtt - high familiarity might produce high effort in first attempt, but in the following attempts there will be less needed
so
- Familiarity
- (1 + CommAtt | Familiarity) or interaction
- Big5
As it’s influencing CommAtt and Effort, it should be there as fixed effect to block the causal path.
We also expect different baselines of effort depending on how, for example, open people are, so that’s why we go for varying intercept.
Do we expect that personality traits affect people differently? Probably not
so
- Big5
- (1 | Big5)
- Expressibility
As it’s influencing CommAtt and Effort, it should be there as fixed effect to block the causal path.
Do we expect different baselines of effort depending on how the concept is expressible? Yes, probably
Do we expect that expressibility could influence effort differently? Not sure
Model 8: full adapted
priors_eff8 <- c(
prior(normal(2.5, 1), class = "Intercept"), # Prior for the intercept (baseline Effort)
prior(normal(0, 1.5), class = "b", coef = "CommAtt2M1"), # Prior for CommAtt level 2
prior(normal(0, 1.5), class = "b", coef = "CommAtt3M2"), # Prior for CommAtt level 3
prior(normal(0, 0.5), class = "b", coef = "Modality1"), # Prior for Modality level 1
prior(normal(0, 0.5), class = "b", coef = "Modality2"), # Prior for Modality level 2
prior(normal(0, 0.5), class = "b", coef = "Familiarity"),
prior(normal(0, 0.5), class = "b", coef = "Expressibility_z"),
prior(normal(0, 0.5), class = "b", coef = "Big5"),
prior(normal(0, 0.5), class = "sd", group = "Participant"), # Random effect for Participant (Intercept)
prior(normal(0, 0.5), class = "sd", group = "Concept"), # Random effect for Concept (Intercept)
prior(normal(0, 0.5), class = "sd", group = "Familiarity"), # Random effect for Familiarity (Intercept)
prior(normal(0, 0.5), class = "sd", group = "Expressibility_z"), # Random effect for Expressibility (Intercept)
prior(normal(0, 0.5), class = "sd", group = "TrialNumber_c"), # Random effect for Trial Number (Intercept)
prior(exponential(1), class = "sigma"), # Prior for residuals
prior(lkj(4), class = "cor")
)
fit_eff8 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 | Big5) + (1 | Expressibility_z) + (1 | TrialNumber_c),
data = final_data,
prior = priors_eff8,
cores = 4)
saveRDS(fit_eff8, here("06_Modeling", "models", "fit_eff8.rds"))
beep(5)fit_eff8 <- readRDS(here("06_Modeling", "models", "fit_eff8.rds"))
# summary
summary(fit_eff8)## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 | Concept) + (1 + CommAtt | Familiarity) + (1 | Big5) + (1 | Expressibility_z) + (1 | TrialNumber_c)
## Data: final_data (Number of observations: 5055)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~Big5 (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.01 0.01 0.00 0.04 1.00 1736 1639
##
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.02 0.01 0.00 0.05 1.00 1140 1639
##
## ~Expressibility_z (Number of levels: 63)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.05 0.02 0.01 0.08 1.01 582 601
##
## ~Familiarity (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.02 0.01 0.00 0.05 1.06 53
## sd(CommAtt2M1) 0.22 0.11 0.00 0.34 1.53 7
## sd(CommAtt3M2) 0.29 0.15 0.01 0.44 1.53 7
## cor(Intercept,CommAtt2M1) 0.27 0.34 -0.45 0.81 1.13 23
## cor(Intercept,CommAtt3M2) -0.25 0.34 -0.79 0.47 1.13 22
## cor(CommAtt2M1,CommAtt3M2) -0.73 0.41 -0.98 0.37 1.53 7
## Tail_ESS
## sd(Intercept) 775
## sd(CommAtt2M1) 33
## sd(CommAtt3M2) 32
## cor(Intercept,CommAtt2M1) 170
## cor(Intercept,CommAtt3M2) 168
## cor(CommAtt2M1,CommAtt3M2) 29
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.02 0.01 0.00 0.05 1.08 38
## sd(CommAtt2M1) 0.09 0.11 0.00 0.32 1.53 7
## sd(CommAtt3M2) 0.12 0.15 0.00 0.42 1.53 7
## cor(Intercept,CommAtt2M1) 0.12 0.35 -0.55 0.77 1.16 18
## cor(Intercept,CommAtt3M2) -0.11 0.36 -0.77 0.56 1.17 18
## cor(CommAtt2M1,CommAtt3M2) -0.29 0.47 -0.97 0.53 1.53 7
## Tail_ESS
## sd(Intercept) 73
## sd(CommAtt2M1) 34
## sd(CommAtt3M2) 30
## cor(Intercept,CommAtt2M1) 54
## cor(Intercept,CommAtt3M2) 58
## cor(CommAtt2M1,CommAtt3M2) 29
##
## ~TrialNumber_c (Number of levels: 54)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.03 0.01 0.00 0.06 1.00 1223 1532
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.72 0.01 2.69 2.75 1.00 3194 2949
## CommAtt2M1 1.60 0.03 1.54 1.67 1.00 1625 2534
## CommAtt3M2 -2.13 0.04 -2.21 -2.04 1.00 1843 2688
## Modality1 -0.55 0.03 -0.62 -0.48 1.00 2696 2881
## Modality2 0.30 0.03 0.24 0.36 1.00 2952 3221
## Big5 1.15 0.04 1.07 1.22 1.01 833 1880
## Familiarity 1.17 0.04 1.09 1.25 1.01 758 1738
## Expressibility_z 0.39 0.01 0.37 0.41 1.00 3393 2397
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.62 0.01 0.61 0.64 1.00 5500 2915
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# not that there are troubles with convergence! (default priors)
plot(fit_eff8)
# looks scheisse, divergent chains
plot(conditional_effects(fit_eff8), points = TRUE)
# nicht so gut
pp_check(fit_eff8, type = "dens_overlay")## Using 10 posterior draws for ppc type 'dens_overlay' by default.
# stays the same
pp_check(fit_eff8, type = "error_scatter_avg")## Using all posterior draws for ppc type 'error_scatter_avg' by default.
# Let’s do the same model, but familiarity will not have varying slope
Model 9: full adapted 2
priors_eff9 <- c(
# Prior for the intercept
set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
# Priors for the main effects using Student's t-distribution
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z"),
set_prior("gamma(2,1)", class = "sd", group = "TrialNumber_c"),
set_prior("gamma(2, 1)", class = "sd", group = "Participant"), # let's try gamma because it limits the negative values a bit better
set_prior("gamma(2, 1)", class = "sd", group = "Concept"),
set_prior("gamma(2, 1)", class = "sd"),
# Prior for residual standard deviation (sigma)
#set_prior("normal(0, 2.5)", class = "sigma", lb = 0),
set_prior("gamma(2, 1)", class = "sigma"),
prior("lkj(2)", class = "cor") # lkj assumes no extreme correlation
)
fit_eff9 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 | Concept) + (1 | Expressibility_z) + (1 | TrialNumber_c),
data = final_data,
prior = priors_eff9,
cores = 4)
saveRDS(fit_eff9, here("06_Modeling", "models", "fit_eff9.rds"))
beep(5)Ok now I am making a bit better sense of it
so we definitely have our main predictor CommAtt
We have other predictors that we believe moderate the effect - that is modality, familiarity, expressibility, big5
we want to have also varying intercept and slope for some
participants - because we believe that their baseline effort varies, and that the effect of commatt on eff might vary across them
same for concepts
we want to have varying intercept for TrialNumber because there might be variation between earlier and later performances (either learning, or fatigue)
Making sense of priors
We need to look a bit closer at our priors so that it’s reasonable
priors_eff9 <- c(
set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z"),
set_prior("gamma(2,1)", class = "sd", group = "TrialNumber_c"),
set_prior("gamma(2, 1)", class = "sd", group = "Participant"),
set_prior("gamma(2, 1)", class = "sd", group = "Concept"),
set_prior("gamma(2, 1)", class = "sd"),
set_prior("gamma(2, 1)", class = "sigma"),
set_prior("lkj(2)", class = "cor") # lkj assumes no extreme correlation
)
fit_eff9_priors <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c),
data = final_data,
prior=priors_eff9,
family = gaussian,#
sample_prior = 'only',
cores = 4)## Compiling Stan program...## Start samplingpp_check(fit_eff9_priors, type = "dens_overlay", ndraws = 100)
pp_check(fit_eff9_priors,
type = "stat",
stat = "mean",
prefix = "ppd") +
coord_cartesian(xlim = c(-20, 20)) +
ggtitle("Prior predictive distribution of means")## Using all posterior draws for ppc type 'stat' by default.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
pp_check(fit_eff9_priors,
type = "stat",
stat = "min",
prefix = "ppd") +
coord_cartesian(xlim = c(-50, 20)) +
ggtitle("Prior predictive distribution of means")## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# this is still of, the minimum value should be 0
pp_check(fit_eff9_priors,
type = "stat",
stat = "max",
prefix = "ppd") +
coord_cartesian(xlim = c(-20, 50)) +
ggtitle("Prior predictive distribution of means")## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Let’s try priors with exponential(1) for sd, whether it will fix the negative values
priors_eff_exp <- c(
set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z"),
set_prior("exponential(2)", class = "sd", group = "TrialNumber_c"),
set_prior("exponential(2)", class = "sd", group = "Participant"),
set_prior("exponential(2)", class = "sd", group = "Concept"),
set_prior("exponential(2)", class = "sd"),
set_prior("exponential(1)", class = "sigma"),
set_prior("lkj(2)", class = "cor") # lkj assumes no extreme correlation
)
fit_eff9_priors_exp <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c),
data = final_data,
prior=priors_eff_exp,
family = gaussian,#
sample_prior = 'only',
cores = 4)## Compiling Stan program...## Start samplingpp_check(fit_eff9_priors_exp, type = "dens_overlay", ndraws = 100)
pp_check(fit_eff9_priors_exp,
type = "stat",
stat = "mean",
bins = 50,
prefix = "ppd") +
coord_cartesian(xlim = c(-20, 20)) +
ggtitle("Prior predictive distribution of means")## Using all posterior draws for ppc type 'stat' by default.
pp_check(fit_eff9_priors_exp,
type = "stat",
stat = "min",
prefix = "ppd") +
coord_cartesian(xlim = c(-50, 20)) +
ggtitle("Prior predictive distribution of minimal values")## Using all posterior draws for ppc type 'stat' by default.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# this is still of, the minimum value should be 0
pp_check(fit_eff9_priors_exp,
type = "stat",
stat = "max",
prefix = "ppd") +
coord_cartesian(xlim = c(-20, 50)) +
ggtitle("Prior predictive distribution of maximal values")## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Let’s try to tighten the priors even more , especially the sds that are
still pushing the posterior to negative values (which is not
possible)
priors_eff_t <- c(
set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z"),
set_prior("normal(0.5,0.1)", class = "sd", group = "TrialNumber_c"),
set_prior("normal(0.5,0.1)", class = "sd", group = "Participant"),
set_prior("normal(0.5,0.1)", class = "sd", group = "Concept"),
set_prior("normal(1,0.1)", class = "sd"),
set_prior("normal(0.5,0.25)", class = "sigma"),
set_prior("lkj(2)", class = "cor") # lkj assumes no extreme correlation
)
fit_eff9_priors_t <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c),
data = final_data,
prior=priors_eff_t,
family = gaussian,#
sample_prior = 'only',
cores = 4)## Compiling Stan program...## Start samplingpp_check(fit_eff9_priors_t, type = "dens_overlay", ndraws = 100)
pp_check(fit_eff9_priors_t,
type = "stat",
stat = "mean",
bins = 50,
prefix = "ppd") +
coord_cartesian(xlim = c(-10, 10)) +
ggtitle("Prior predictive distribution of means")## Using all posterior draws for ppc type 'stat' by default.
# this look okay
pp_check(fit_eff9_priors_t,
type = "stat",
stat = "min",
prefix = "ppd") +
coord_cartesian(xlim = c(-15, 10)) +
ggtitle("Prior predictive distribution of minimal values")## Using all posterior draws for ppc type 'stat' by default.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# this is better but we still see some 0s
pp_check(fit_eff9_priors_t,
type = "stat",
stat = "max",
prefix = "ppd") +
coord_cartesian(xlim = c(-10, 15)) +
ggtitle("Prior predictive distribution of maximal values")## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# this looks reasonableCheck out: https://paulbuerkner.com/brms/reference/set_prior.htmlNow let’s exchange the priors for the tighter ones
Model 10: full with mildly informative priors
fit_eff10 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c),
data = final_data,
prior=priors_eff_t,
family = gaussian,
iter = 4000,
cores = 4)
saveRDS(fit_eff10, here("06_Modeling", "models", "fit_eff10.rds"))
# warnings with ESS for correlation coefficients!fit_eff10 <- readRDS(here("06_Modeling", "models", "fit_eff10.rds"))
# SUMMARY
summary(fit_eff10)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c)
## Data: final_data (Number of observations: 5055)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.04 0.02 0.01 0.09 1.00 2037
## sd(CommAtt2M1) 0.28 0.06 0.17 0.42 1.00 1824
## sd(CommAtt3M2) 0.26 0.07 0.15 0.41 1.00 2139
## cor(Intercept,CommAtt2M1) 0.25 0.37 -0.53 0.85 1.01 625
## cor(Intercept,CommAtt3M2) -0.20 0.38 -0.84 0.58 1.00 676
## cor(CommAtt2M1,CommAtt3M2) -0.92 0.07 -0.99 -0.71 1.00 3026
## Tail_ESS
## sd(Intercept) 1928
## sd(CommAtt2M1) 3836
## sd(CommAtt3M2) 3584
## cor(Intercept,CommAtt2M1) 1090
## cor(Intercept,CommAtt3M2) 1397
## cor(CommAtt2M1,CommAtt3M2) 4946
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.06 0.02 0.03 0.09 1.01 1067
## sd(CommAtt2M1) 0.32 0.03 0.27 0.38 1.00 3558
## sd(CommAtt3M2) 0.41 0.04 0.34 0.49 1.00 2883
## cor(Intercept,CommAtt2M1) 0.82 0.13 0.47 0.97 1.01 467
## cor(Intercept,CommAtt3M2) -0.82 0.13 -0.97 -0.46 1.02 468
## cor(CommAtt2M1,CommAtt3M2) -0.97 0.02 -1.00 -0.93 1.00 4931
## Tail_ESS
## sd(Intercept) 914
## sd(CommAtt2M1) 5393
## sd(CommAtt3M2) 4748
## cor(Intercept,CommAtt2M1) 357
## cor(Intercept,CommAtt3M2) 323
## cor(CommAtt2M1,CommAtt3M2) 6177
##
## ~TrialNumber_c (Number of levels: 54)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.04 0.01 0.01 0.07 1.00 2168 2276
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.72 0.02 2.69 2.75 1.00 2992 3779
## CommAtt2M1 1.53 0.07 1.38 1.67 1.00 1816 3049
## CommAtt3M2 -2.05 0.08 -2.19 -1.89 1.00 2016 2732
## Modality1 -0.54 0.03 -0.60 -0.49 1.00 6604 6447
## Modality2 0.30 0.03 0.25 0.35 1.00 8200 6517
## Big5 1.08 0.04 0.99 1.16 1.01 1662 2427
## Familiarity 1.09 0.05 1.00 1.18 1.01 1284 1397
## Expressibility_z 0.39 0.01 0.37 0.41 1.00 4057 5432
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.62 0.01 0.61 0.63 1.00 11144 5631
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# PLOT COEFFS
plot(fit_eff10)
plot(conditional_effects(fit_eff10), points = TRUE)
# PP CHECK
pp_check(fit_eff10, type = "dens_overlay", ndraws = 100)
# the fit looks quite ok
## PP CHECK WITH SUMMARY STATS
pp_check(fit_eff10,
type = "stat",
stat = "mean",
bins = 50)## Using all posterior draws for ppc type 'stat' by default.
pp_check(fit_eff10,
type = "stat",
stat = "min",
bins = 50) # ok here we see it's looking for values where it should not## Using all posterior draws for ppc type 'stat' by default.
pp_check(fit_eff10,
type = "stat",
stat = "max",
bins = 50) # here it's also quite off## Using all posterior draws for ppc type 'stat' by default.
## PP CHECK FOR RESIDUALS
pp_check(fit_eff10, type = "error_scatter_avg")## Using all posterior draws for ppc type 'error_scatter_avg' by default.
Where to next
What are potential interactions?
- Modality*Expressibility (only highly expressive concepts might give modality the power to affect effort)
- Big5*Modality (only in vocal modality traits could matter)
- Big5 and/or Familiarity * CommAtt
Nesting
Dyad:Participant Modality:Concept Modality:Expressibility -
possibly change family to student , family = student,